Investigating Intrinsically Disordered Proteins With Brownian Dynamics

Intrinsically disordered proteins (IDPs) have recently become systems of great interest due to their involvement in modulating many biological processes and their aggregation being implicated in many diseases. Since IDPs do not have a stable, folded structure, however, they cannot be easily studied with experimental techniques. Hence, conducting a computational study of these systems can be helpful and be complementary with experimental work to elucidate their mechanisms. Thus, we have implemented the coarse-grained force field for proteins (COFFDROP) in Browndye 2.0 to study IDPs using Brownian dynamics (BD) simulations, which are often used to study large-scale motions with longer time scales and diffusion-limited molecular associations. Specifically, we have checked our COFFDROP implementation with eight naturally occurring IDPs and have investigated five (Glu-Lys)25 IDP sequence variants. From measuring the hydrodynamic radii of eight naturally occurring IDPs, we found the ideal scaling factor of 0.786 for non-bonded interactions. We have also measured the entanglement indices (average C α distances to the other chain) between two (Glu-Lys)25 IDP sequence variants, a property related to molecular association. We found that entanglement indices decrease for all possible pairs at excess salt concentration, which is consistent with long-range interactions of these IDP sequence variants getting weaker at increasing salt concentration.

The system of equations can be written as where r is the vector of all particle positions after a time step without any constraints, Φ is the equation describing the constraints, M is the collective mobility matrix of the particles, J is a generalized impulse, with one component for each constraint, and the matrix C converts the generalized impulses into impulses (force times time step) on the particles. Because corrective impulses should be applied to the particles in the direction by which they most violate the constraint, the matrix C is defined as Because J has units of impulse and M has units of velocity over force, the matrix C must be dimensionless and the constraints Φ must have units of length.
The system of equations, Eq. S1 can be solved by Newton-Raphson (NR) iteration by repeatedly solving for F and updating the positions by r ← r + M · C · J (S4) In the simulations here, the bond length constraint is used: where subscripts i(1) and i(2) refer to the two particles of the i th constraint, and R i is the desired distance. One can also devise other constraints, such as coplanarity and collinearity. The iterations are performed until each constraints is satisfied within tolerance ϵ: In the implementation, the matrix in Eq. S3 is represented as a skyline matrix, which is a banded matrix with the bandwidth varying with the row. It is solved using a skyline Cholesky solver, which does not introduce any new non-zero elements into the matrix during factorization. The implementation automatically computes the bandwidth for each row based on the connectivity of the constraints with the particles, and the bandwidth of the matrix is minimized by numbering the constraints in order of their place along the chain.
In practice, this NR method often does not converge unless the constraints are already close to being satisfied (in which case it converges very quickly). Two more layers are added to the solver in order to make it robust and stable.
First, if a NR step increases the root-mean-square of the constraint violations from Eq. S6, or the Cholesky factorization fails, then an iteration of a SHAKE-like method (Ryckaert et al., 1977) for each constraint is performed where a i is the radius of bead i and α is a relaxation factor that starts as 1. If the iteration fails to decrease the RMS violation, then the bead positions are reset and α is halved. Once the RMS violation has decreased, the NR iterations continue. The NR iterations are fast but uncertain, while the SHAKE iterations are slow but sure.
Before the NR-SHAKE iterations are performed, another layer is used to ensure stability. If any of the constraints for the new positions from the BD step are violated by more than a factor of 0.1, then the bead positions are reset, and pushed along the directions from the initial bead positions to the final bead positions from the unconstrained step.
where r bi is the beginning position, r f i is the final position, and r ic = r bi at the start. A value of parameter λ is found, using bisection, for which the positions {r i } give a maximum violation of 0.1. The constraints are solved, and the positions r ci are set to the newly constrained positions. This procedure is repeated until all the values of λ from the steps add up to 1.

APPENDIX B: HYDRODYNAMICS ALGORITHM
Following work by Elcock (Elcock, 2013), we have implemented an approximation which splits the hydrodynamic interactions into local and detailed effects (in this case, intramolecular), and distant and coarser effects (intermolecular). We have also included the stochastic forces in a way that preserves the fluctuation-dissipation theorem. The required time of the simplest algorithm, which computes the Cholesky decomposition of the Rotne-Prager-Yamakawa (RPY) mobility tensor, scales as O(n 3 ) where n is the number of spherical hydrodynamic "beads". Algorithms with better scaling do exist, but the more modestly sized systems studied here are probably not large enough to benefit from the advanced algorithms.
Suppose we have m molecules, or separate units, and the i th molecule has n i beads and a 3n i × 3n i mobility matrix M ii . We replace the mobility matrix M ij between separate molecules with a coarse-grained treatment. The algorithm begins with separately computing the Cholesky decomposition for each M ii . For example, for two equally-sized molecules, two separate decompositions, neglecting the intermolecular terms, will be about four times faster than the decomposition of the full matrix. Next, the 3 × 3 mobility matrix D ci of each molecule as a rigid body, and its mobility center c i , is computed by (insert here). From the mobility of each molecule, an effective radius is computed: where µ is the solvent viscosity, and a 3m × 3m mobility matrix M CC is computed from the RPY formula using the radii {a i } and the distances between the centers {c i }.
The drift velocity of each molecule, as a rigid body, is computed using where U 1 is the 3 × 3 identity matrix and U p is U 1 stacked upon itself p times. The forces and velocities of the beads of Molecule i are the 3n i vectors F i and v. This is equivalent to multiplying the sum of the forces on each molecule by its mobility matrix M Ci , and assigning the resulting velocity to each bead.
Next, local drift velocity components of each molecule are computed and added to the velocity components from Eq. S13. The p × p identity matrix is I p . This is equivalent to subtracting the average force on each bead, multiplying the local mobility matrix M i , and subtracting the average velocity from each bead. The complete computation of the drift velocity can be summarized by where the actions of matrices M C and M F can be summarized by Eqs. S13 and S14, respectively. Both matrices are symmetric and positive definite, and their square roots can be computed as where L CC and L ii are the Cholesky decompositions of the matrices M CC and M ii . It can be shown, that for stochastic vectors w C , w F , and w of independent Gaussian distributions of zero mean and unit variance, the stochastic quantity has the same properties (such as covariance) as the quantity which is required to satisfy the fluctuation-dissipation theorem. Thus, the stochastic term of the Brownian dynamics can be computed using Eq. S18. As discussed above, for several molecules or subunits, computing several smaller Cholesky decompositions {L ii } will be quicker, and computing the decomposition of M CC should not take much time at all, provided the number of beads is much greater than the number of molecules. Future work will include the addition of torques and rotational motion to the formulation for greater accuracy, since extensions to the RPY tensor now exist for rotational motion (Zuk et al., 2014).   Figure S1. Autocorrelation functions of end-to-end C α distance, end-to-middle C α distance, and middleto-end C α distance for (Glu-Lys) 25 sv=10 at NaCl 15 mM (normal concentration) with scaling factor 1.0. . Autocorrelation functions of end-to-end C α distance, end-to-middle C α distance, and middleto-end C α distance for (Glu-Lys) 25 sv=10 at NaCl 125 mM (excess salt concentration) with scaling factor 1.0. . Autocorrelation functions of end-to-end C α distance, end-to-middle C α distance, and middleto-end C α distance for (Glu-Lys) 25 sv=10 at NaCl 15 mM (normal concentration) with scaling factor 0.825.  Figure S4. Autocorrelation functions of end-to-end C α distance, end-to-middle C α distance, and middleto-end C α distance for (Glu-Lys) 25 sv=10 at NaCl 125 mM (excess salt concentration) with scaling factor 0.825.  Figure S8. Autocorrelation functions of end-to-end C α distance, end-to-middle C α distance, and middleto-end C α distance for (Glu-Lys) 25 sv=15 at NaCl 125 mM (excess salt concentration) with scaling factor 0.825.  Figure S9. Autocorrelation functions of end-to-end C α distance, end-to-middle C α distance, and middleto-end C α distance for (Glu-Lys) 25 sv=20 at NaCl 15 mM (normal concentration) with scaling factor 1.0.